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The method of iterated conformal maps is developed for quasi-static fracture of brittle materials, 
for all modes of fracture. Previous theory, that was relevant for mode III only, is extended here to 
mode I and II. The latter require solution of the bi-Laplace rather than the Laplace equation. For 
all cases we can consider quenched randomness in the brittle material itself, as well as randomness 
in the succession of fracture events. While mode III calls for the advance (in time) of one analytic 
function, mode I and II call for the advance of two analytic functions. This fundamental difference 
creates different stress distribution around the cracks. As a result the geometric characteristics of 
the cracks differ, putting mode III in a different class compared to modes I and II. 



I. INTRODUCTION 



at infinity. These are 



The theory of quasi-static fractures in brittle media 
10, ^ m ^, 1^ calls for solving different equations depend- 
ing on the mode of fracture. In this paper we present an 
approach based on iterated conformal maps which can 
be adapted to solve all three modes of fracture (known 
as mode I, II and III), including the effects of inhomo- 
geneities and randomness of the brittle material itself. 

Basically, the theory of fracture in brittle continuous 
media is based on the equation of motion for an isotropic 
elastic body in the continuum limit Q 



(A-l-^)V(V-u)-f /iV^u 



(1) 



Here u is the field describing the displacement of each 
mass point from its location in an unstrained body and 
p is the density. The constants fi and A are the Lame 
constants. In terms of the displacement field the elastic 
strain tensor is defined as 



dui 



dxi 



(2) 



For the development of a crack the important object is 
the stress tensor, which in linear elasticity is written as 



\5r 



k 



(3) 



When the stress component which is transverse to the 
interface of a crack exceeds a threshold value ctc, the 
crack can develop. When the external load is such that 
the transverse stress exceeds only slightly the threshold 
value, the crack develops slowly, and one can neglect the 
second time-derivative in Eq. (Q). This is the quasi-static 
limit, in which after each growth event one needs to re- 
calculate the strain field by solving the Lame equation 



(A -1- Ai)V(V • u) + ^V^u = 



(4) 



The three "pure" modes of fracture that can be consid- 
ered are determined by the boundary conditions, or load. 



o-i2;(oo) = ;a-yj^(oo) = (Too ;o-a;y(oo) = Model (5) 
(Txxioo) = ; cTyy (oo) = ;crxy{oo) = tJoo Mode II (6) 

We will study the fracture patterns of these two modes 
in 2-dimensional materials. Mode III calls for a third 
dimension z, since 



o'zyiy ±oo) — Goo Mode III. 



(7) 



Such an applied stress creates a displacement field 
Uz{x,y), Ux = 0, Uy = in the medium. Thus, in spite 
of the third dimension, the calculation of the strain and 
stress tensors remain two dimensional. Nevertheless, the 
equations to be solved in mode III and modes I and II 
are different. In mode III fracture V • u = 0, and the 
Lame equation reduces to Laplace's equation 



(8) 



and therefore is the real part. Re x(2;), of an analytic 
function x(^); 



Xiz) = Uz{x,y) +i^z{x,y) , 



(9) 



where z = x + iy. The boundary conditions far from the 
crack and on the crack interface can be used to find this 
analytic function. On the other hand, for mode I and 
mode II fractures in plane elasticity one introduces 0| 
the Airy potential U{x, y) such that 



9y2 '--y Q^Qy '^''yy - q^2 ■ (10) 

The Airy potential U solves the bi-Laplacian equation ||] 

AAf7(a;, y) = . (11) 

The solution of the bi-Laplacian equation can be written 
in terms of two analytic functions <f)(^z) and r\{z) as 



U{x,y) = Re[zv3(z) + 77(z)] 



(12) 
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This difference requires therefore a separate discussion of 
mode III and modes I and II. 

The problem of quasi-static crack propagation is diffi- 
cuh not only because it is hard to solve Eq. (^) for an 
arbitrarily shaped crack. Another source of difficulty is 
that the equation does not dictate how to propagate a 
crack when the stress tensor exceeds the threshold value 
CTc. In this paper we consider only 2-dimensional, or ef- 
fectively 2-dimensional (i.e thin slabs) brittle materials in 
a;, y. We can then describe a crack of arbitrary shape by 
its interface x{s), where s is the arc length which is used 
to parameterize the contour. We will use the notation 
(t, n) to describe respectively the transverse and normal 
directions at any point on the two-dimensional crack in- 
terface. The literature is quite in agreement that the 
velocity of propagation of the crack has a normal com- 
ponent which is some function of (Ttt{s) — (Jc for mode I 
and II, and of |cr2t(s)| — CTc for mode III. In both cases CTc 
is a measure of the strength of the material, and fracture 
occurs only if the local stress tensor at the boundary of 
the crack exceeds this quantity (which can also be a ran- 
dom function of position). There is hardly a consensus 
however on what that function is. The simplest choice 
1^, is a linear function, 

Vn{s) = aAa = a{cru{s) - (Jc{s)) , Mode I,II ,(13) 
Vn{s) = aAa = a{\(j^t{s)\ - (Jc{s)) , mode III (14) 

when Act > 0, and Vn{s) = otherwise. Other velocity 
laws are possible |^ . In our study of mode III fracture we 
will examine also a quadratic and an exponential velocity 
law: 



„Q(k^t(s)|-o-c(s)) 



mode III 



mode III 



(15) 
(16) 



It is important to study these variants of the velocity law 
to ascertain the degree of universality of the geometric 
characteristics of the resulting cracks. One of our results 
is that these characteristics may depend on the velocity 
law. While this may be a disappointment from the point 
of view of fundamental physics, it may help to identify 
the correct physical mechanisms of fractures in different 
media. The lack of universality is even more obvious 
when we add quenched noise, or random values of ads). 
The geometric characteristics of the cracks may depend 
on the probability distribution of random values of ads). 
Again this may give a handle on the characterization of 
inhomogeneous brittle materials. 

At any point in time there can be more than one po- 
sition s on the interface for which w„(s) does not van- 
ish. We choose the next growth position randomly with 
a probability proportional to f„(s) [0, ||]. There we ex- 
tend the crack by a fixed area of the size of the "process 
zone" (and see below for details). This is similar to Dif- 
fusion Limited Aggregation (DLA) in which a particle is 
grown with a probability proportional to the gradient of 
the field. One should note that another model could be 



derived in which all eligible fracture sites are grown si- 
multaneously, growing a whole layer whose local width 
is Vn{s). This would be more akin to Laplacian growth 
algorithms, which in general give rise to clusters in a dif- 
ferent universality class than DLA jl^, |ll| . 

In Sect. II we discuss the growth al gori thm in terms 



III] this method 
A prelim- 



of iterated confornial maps. In Sect 
is applied to mode III quasi-static fracture 
inary report of the method for this case was presented 
in In Sect. IV we present new results including 

the consequences of the different velocity laws (|l^) and 
(p^), and those of quenched randomness. We discuss the 
geometric properties of the fracture patterns, including 
issues of roughening and exponents. We point out that 
the roughening exponents are not always well defined, 
since the fracture patterns do not have stationary ge- 
ometric characteristics. There is an increased tendency 
for ramification as the fracture develops. This is reflected 
in an apparent increase in the roughening exponents of 
the backbone of the pattern. In Sect. Mwe discuss the 
theory of mode I and II fracture. Sect. VI presents the 
results. We will see that the fracture patterns in mode I 
and II are much less rough then in mode III (for the same 
velocity law), in agreement with the a nalysis of |l^. We 
will conclude the paper in Sect. VIL The main conclu- 



sion is that mode III results in cracks whose geometric 
characteristics are in a different class than modes I and 
II. The former creates cracks that exhibit a cross over 
in the roughening exponent from about 0.5 to a higher 
scaling exponent on the larger scales. In contrast, modes 
I and II create cracks that are not rough on the large 
scales. Quenched randomness may affect the geometry 
of the cracks as is exemplified and discussed in this pa- 
per. 



II. THE METHOD OF ITERATED 
CONFORMAL MAPS FOR FRACTURE 

The direct determination of the strain tensor for an 
arbitrary shaped (and evolving) crack is difficult. We 
therefore proceed by turning to a mathematical complex 
plane uj, in which the crack is forever circular and of unit 
radius. Next invoke a conformal map z = $(")(w) that 
maps the exterior of the unit circle in the mathemati- 
cal plane uj to the exterior of the crack in the physical 
plane z, after n growth steps. The conformal map will be 
univalent by construction, and we can write its Laurent 
expansion in the form 



$(")(,, 



(") I 







f^1>/uj + f[":^/ 



(17) 



For all modes of fracture we take <I>*^°^(a;) — uj, and the 
iterative dynamics calls for the calculation of the trans- 
verse component of the stress tensor on the boundary of 
the crack. The arclength position s in the physical do- 
main is mapped by the inverse of onto a position the 
unit circle uj = exp{i9). We will be able to compute the 
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stress tensor on the boundary of the crack in the phys- 
ical domain by performing the calculation on the unit 
circle. In other words we will compute att{6) or (Tzt{d) 
on the unit circle in the mathematical plane. The actual 
calculation of this component of the stress tensor differs 
in modes I, II and mode III. We perform the calculation 
iteratively, taking the stress as known for the crack after 
n — 1 fracture events. 

In order to implement the nth cracking event according 
to one of the required velocity laws (|l^-(|l^), we should 
choose potential positions on the interface more often 
when Aa{9) is larger. Consider for example the linear 
velocity law (|l3|). We construct a probability density 
P{9) on the unit circle e** which satisfies 



P{9) 



|$'(n-i)(e'9)|A(T(6>)9(Ag(6l)) 
/o'"|$'("-i)(e»«)|Aa(^)e(A(T(0))d^ 



(18) 



where Q{Aa{9)) is the Heaviside function, and 
1$ ("~^)(e'^)| is simply the Jacobian of the transforma- 
tion from mathematical to physical plane. The next 
growth position, 9n in the mathematical plane, is cho- 
sen randomly with respect to the probability P{9)d9. At 
the chosen position on the crack, i.e. z = <I)("~i)(e'^"), 
we want to advance the crack with a region whose area is 
the typical process zone for the material that we analyze. 
According to ||] the typical scale of the process zone is 
K^/(j1, where K \s a. characteristic fracture toughness 
parameter. Denoting the typical area of the process zone 
by Ao, we achieve growth with an auxiliary conformal 
map (f>\^_g^{uj) that maps the unit circle to a unit circle 
with a bump of area A„ centered at e**" . An example of 
such a map is given by p^ : 



2w 



1 2 1-A 
l + w + w\ l + ^ - 



wl + X 



1/2 



- 1 



(19) 



(20) 



Here the bump has an aspect ratio a, < a < 1. In our 
work below we use a = 1/2. To ensure a fixed size step 
in the physical domain we choose 



An 



An 



(21) 



|$(«-i)'(gie„)|2 

Finally the updated conformal map <I''^") is obtained as 
$(")(c.) = $("-i)(0,^_,Jc^)) . (22) 

The recursive dynamics can be represented as itera- 
tions of the map </)A„,e„(w), 



$(")(w) = 0Ai,ei ° (1^X2,92 O ...O (f)x^M„{^^) 



(23) 



Every given fracture pattern is determined completely by 
the random itinerary 




FIG. 1: Upper panel: a typical mode III fracture pattern 
that is obtained from iterated conformal maps. What is seen 
is the boundary of the fractured zone, which is the mapping of 
the unit circle in the mathematical domain onto the physical 
domain. Notice that the pattern becomes more and more 
ramified as the the fracture pattern develops. This is due to 
the enhancement of the stress field at the tips of the growing 
pattern. Lower panel: the backbone of the fracture pattern. 
This is the projection onto the x-y plane of the experimentally 
observed boundary between the two parts of the material that 
separate when the fracture pattern hits the lateral boundaries. 



III. MODE III QUASI-STATIC FRACTURE 

In this section we discuss how to compute the stress 
tensor when the load is mode HI, using the method of 
iterated conformal maps. The first step is the deter- 
mination of the boundary conditions that the analytic 
function (^ needs to satisfy. 



A. Boundary conditions in mode III 

Far from the crack as y ^ ±cx) we know cr^j, —^ (Too or 
using the stress-strain relationships Eq. (||) we find that 
Uz ~ [aoo/ ^J■]y■ Thus the analytic function must have the 
form 

x{z) -i[aoo/^\z as \z\ ^ oo . (24) 
Now on the boundary of the crack the normal stress 
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vanishes, i.e. 



= (Tzn(s) = dnUz = ~dt£,z 



(25) 



This means that is constant on the boundary. We 
choose the gauge = 0, which in turn is a boundary 
condition making the analytic function xC^) real on the 
boundary of the crack: 



(26) 




lOCX) 



B. The stress tensor for mode III 

Following the basic strategy we consider now a circular 
crack in the mathematical domain. The strain field for 
such a crack is well known Q , being the real part of the 
fimction x^^"^ (^) where 



FIG. 2: h{r) averaged over all the backbone and over 20 frac- 
ture patterns each of which of 10 000 fracture events. There 
is a cross-over between a scaling law with roughness exponent 
0.49 ± 0.08 to an exponent of 0.70 ± 0.05 



IV. RESULTS FOR MODE III 



X^°H^) = -^Woo/fAi^-l/u;) 



(27) 



This is the unique analytic function obeying the bound- 
ary conditions x'^H^) ~^ — ^[coo/m]^ as |ci;| — > oo, while 
on the unit circle x*^") (exp i6') = x^'^H^xpi^)*- To find 
the corresponding function in the physical plane is par- 
ticularly easy for mode 111. Since the real part of the 
fimction xi^) is analytic, it satisfied Laplace's equation 
automatically. We only need to make sure that it sat- 
isfies the boundary conditions. However, if we have a 
good solution in the mathematical plane, we need just to 
compose it with an analytic function that takes us from 
the physical to the mathematical plane. The required 
analytic function x'"^^) given by the expression 

X(")(z) = -z[i^i^"Vo,^]($(")"\z) - l/<D(")"\z)) . 

(28) 

(If $(") is conformal, is analytic by definition). 

From this we should compute now the transverse stress 
tensor: 



o'zt(s) = fi dtUz = 



ax(")($(")(e''')) 39 



|$'(")(e^«)| 

(n) cos 9 



|$'(")(e*«)| 



(29) 



on the boundary. Eqs.(p9|) together with ( p3| ) offer an 
analytic expression for the transverse stress field at any 
stage of the crack propagation. 



A. Linear velocity law 

Fig. |l| exhibits in the upper panel a typical fracture 
pattern that is obtained with this theory, with aoo = 1, 
after 10 000 growth events. The threshold value of ctc for 
the occurrence of the first event (cf. Eq.(|29|) is (Tc = 2. 
We always implement the first event. For the next growth 
event the threshold is dc = 2.34315.... We thus display 
in Fig.l a cluster obtained with CTc = 2.00, to be close 
to the quasi-static limit. Note that here we could opt 
to represent a disorde red m aterial by a random value of 
(Tc, and see Subsect. |IV^ . With fixed Ctc, one should 
observe that as the pattern develops, the stress at the 
active zone increases, and we get progressively away from 
the quasi-static limit. Indeed, as a result of this, for fixed 
boundary conditions at infinity, there are more and more 
values of 9 for which Eq.(|8|) does not prohibit growth. 

Since the tips of the patterns are mapped by to 
larger and larger arcs on the unit circle, the support of 
the probability P{9) increases, and the fracture pattern 
becomes more and more ramified as the process advances. 
The geometric characteristics of the fracture pattern are 
not invariant to the growth. For this reason it makes little 
sense to measure the fractal dimension of the pattern; 
this is not a stable characteristic, and it will change with 
the growth. On the other hand, we should realize that 
the fracture pattern is not what is observed in typical 
experiments. When the fracture hits the boundaries of 
the sample, and the sample breaks into two parts, all 
the side-branches of the pattern remain hidden in the 
damaged material, and only the backbone of the fracture 
pattern appears as the surface of the broken parts. The 
backbone does not suffer from the geometric variability 
discussed above. In the lower panel of Fig. ^ we show the 
backbone of the pattern displayed in the upper panel. 

This backbone is representative of all the fracture pat- 
terns with the linear velocity law. We should note that 
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in our theory there are no lateral boundaries, and the 
backbone shown does not suffer from finite size effects 
which may very well exist in experimental realizations. 

In determining the roughness exponent of the back- 
bone, we should note that a close examination of it re- 
veals that it is not a graph. There are overhangs in 
this backbone, and since we deal with mode III frac- 
turing, the two pieces of material can separate leaving 
these overhangs intact. Accordingly, one should not ap- 
proach the roughness exponent using correlation func- 
tion techniques; these may introduce serious errors when 
overhangs exist Rather, we should measure, for any 
given r, the quantity 

h{r) = (Max{?/(r')}^<r/<a;+r - Mm{y{r'}x^r'<x+r)x ■ 

(30) 

The roughness exponent C is then obtained from 



(31) 



if this relation holds. To get good statistics we average, 
in addition to all x for the same backbone, over many 
fracture patterns. The result of the analysis is shown in 
Fig. i 

We find that the roughness exponent for the backbone 
exhibits a clear cross-over from about 0.5 for shorter dis- 
tances r to about 0.70 for larger distances. Within the 
error bars these results are in a surprising agreement with 
the numbers quoted experimentally, see for example . 
The short length scale exponent of order 0.5 is also in 
agreement with recent simulational results of a lattice 
model [|8j (which is by definition a short length scale 
solution). Bouchaud Q proposed that the cross-over 
stems from transition between slow and rapid fracture, 
from the "vicinity of the depinning transition" to the 
"moving phase" in her terms. Obviously, in our theory 
we solve the quasi-static equation all along, and there is 
no change of physics. In addition, there is no reason to 
expect the experiment to be a pure mode III, and as we 
will see below modes I and II do not show similar rough- 
ening. Nevertheless, as we observed before, the fracture 
pattern begins with very low ramification when the stress 
field exceeds the threshold value only at few positions on 
the fracture interface. Later it evolves to a much more 
ramified pattern due to the increase of the stress fields 
at the tips of the mature pattern. The scaling properties 
of the backbone reflect this cross- over. We propose that 
this effect is responsible for the cross-over in the rough- 
ening exponent of the backbone. On the other hand, this 
non-stationarity in the geometric characteristics should 
be handled with care, since it may mean that there is no 
definite roughening exponent, as it may depend on where 
the analysis is done, near the center of the fracture pat- 
terns or near the edge. We will return to this delicate 
issue after reviewing the results of other velocity laws. 
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FIG. 3: Upper panel; fracture pattern for mode III frac- 
ture with the quadratic law (|l5|), with 10 000 fracture events. 
Lower panel: the backbone of the pattern 



B. Other velocity laws 

It should be stressed that there is no reason to believe 
that the scaling exponents are invariant to the change 
of the velocity law. In Figs. ^, ^ and || we show the 
fracture patterns and their corresponding backbones for 
the quadratic velocity law ( p^ ) and for two different ex- 
ponential laws (|l^. We find that the quadratic law 
makes little difference with respect to the linear law. The 
roughening plot is similar, and the scaling exponents ap- 
pear the same. The exponential velocity law changes the 
degree of ramification, and therefore calls for a careful 
discussion of the roughening plots. Examine the func- 
tion h{r) for the pattern in Fig. (|) (see Fig. |). While 
the small scale roughening exponent of about 0.5 is re- 
produced, it appears that the large scale exponent is now 
higher, about 0.78. The question to be asked therefore is 
whether the scaling exponent is not invariant to the ve- 
locity law. In our opinion this question is ill-posed since 
the scaling exponent itself depends on where is it mea- 
sured. As we said before, the fracture pattern tends to 
become more ramified as it grows. This is reflected in 
the roughening properties. To make this point clearer, 
we have taken the pattern of Fig. |5| as a test case, and 
computed the apparent scaling exponents for short parts 
of the fracture pattern, limiting the maximal value of r to 
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FIG. 4: Upper panel: fracture pattern for mode III fracture 
with the exponential law (^6|), with a = 0.1, with 10 000 
fracture events. Lower panel: the backbone of the pattern 
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FIG. 5: Fracture pattern for mode III fracture with the 
exponential law (^6|), with a — 1, with 10 000 fracture events. 
In this case the fracture pattern and the backbone are the 
same. 



2000. Doing so, we can concentrate on a region near the 
center of the pattern, and on a region near the edge. The 
results of this exercise are presented in Fig. ^ What is 
found is that the apparent scaling exponent depends on 
the region of measurements. Near the center, where the 
pattern is less ramified, the exponent is smaller than near 
the edge where the pattern is more ramified. The average 




FIG. 6: h{r) averaged over 20 fracture patterns with the 
exponential velocity law with a = 1. Each of the patterns 
consists of 10 000 fracture events. There is a cross-over be- 
tween a scaling law with roughness exponent of about 0.50 at 
short length scales to an apparent scaling exponent of about 
0.78. 
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FIG. 7: h{r) averaged over 20 fracture patterns with the 
exponential velocity law with a = 1. In this calculation we 
concentrate on parts of the pattern shown in Fig. ^ one near 
the center and the other near the edge, each consisting of 
r = 2000. The apparent exponents differ, being 0.71 at the 
center and 0.85 near the edge. The average behavior with 
exponent 0.78 seen in Fig. M should be therefore interpreted 
with extra care. 



exponent reported in Fig. |^ which is analogous to what 
is reported in experiments, has therefore a limited value. 
It may not be interpreted as a 'true' scaling exponents. 
Its value may well depend on the actual length of the 
pattern that is investigated. 

We are therefore not in a position to claim that the cor- 
respondence in roughening exponents between the linear 
law and experiments indicates anything about universal- 
ity classes. One needs to ascertain very carefully whether 
measured roughening exponents indicate translationally 
invariant scaling properties. It is in particular useful to 
know whether the observed scaling exponents depends on 
the length of the available fracture pattern. 
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FIG. 8: Upper panel: fracture pattern for mode III fracture 
with the linear velocity law and quenched randomness with 
a flat distribution, amax ~ 15, with 10 000 fracture events. 
Lower panel: the function h{r) after averaging over 20 pat- 
terns. The scaling exponents are about 0.4 and 0.65 for the 
smaller and larger scales, respectively. 

C. quenched disorder 

To study the effect of quenched randomness we assign 
a-priori a random value CTc to every point in the mate- 
rial (with resolution Aq). Not having a clear indication 
from the literature how the randomness of inhomoge- 
neous media should be modeled, we opted for two types 
of quenched randomness. The first takes the numerical 
value of (Tc(s) from a flat distribution, < CTc < <^max 
and the second takes a power law form 

P(o-c) oc a^^ , for tJc > o-mm (32) 

For reasonable values of (Tmax the flat distribution did 
not lead to a qualitative change in the fracture patterns. 
In Fig. H we show the pattern and the function h{r) for 
the case amax = 15. The typical cross over that we see in 
systems without quenched disorder remains here, albeit 
with apparently smaller exponents, of about 0.4 and 0.65. 

On the other hand, a power law distribution of 
quenched randomness may lead to very interesting quali- 
tative change in fracture pattern. While high values of /3 
in ( ^ ) are still in qualitative agreement with all previous 
results (see Fig. ^ with [3 = 2), lower values of (3 lead to 



FIG. 9: Upper panel: Fracture pattern for mode III fracture 
with the linear velocity law and quenched randomness with a 
power-law distribution, {3 = 2, Omin = 2, with 10 000 fracture 
events. Lower panel: the function h{r) after averaging over 
20 patterns. The scaling exponent is about 0.65. 
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FIG. 10: Fracture pattern for mode III fracture with the lin- 
ear velocity law and quenched randomness with a power-law 
distribution, (3 = 1.1, Omin = 0.2. 



a new phenomenon. The availability of very high values 
of dc results in effective blocking for the evolution of the 
fracture. The crack develops along continuous (sometime 
curved) lines, and then it suddenly gains sharp turns. 
In Fig. ( |l0| ) we show the typical patterns obtained for 
[3 = 1.1. It is amusing to note that these patterns are 
reminiscent of what is exhibited in a number of experi- 
ments and see for example the pictures in p^. It is not 
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obvious however how to offer quantitative measures for 
comparison. It appears to the present authors that this 
subject of fracture with quenched randomness deserves 
a careful separate study in which experimental and the- 
oretical methods were combined to gain further insights 
on the questions at hand. 



remaining freedoms are now used to choose ipo = and 
tpi real. Then, using the boundary conditions and 
(pf), we find 



^1 



lUr, 



Mode I , 
Mode II . 



(40) 



THEORY FOR MODE I AND II 



B. The Conformal map and its consequences 



In order to compute the stress tensor at the boundary 
of the crack for modes I and II loading, we turn to the 
solution of Eq.(|ll|). 



A. Boundary conditions and removal of freedoms 

The boundary conditions at infinity are given by Eqs. 
® ' ® ■ The conditions on the boundary of the crack are 



o'xn(s) — '^ynis) = On thc boundary 



(33) 



Using Eq. ( pij| ) these boundary condition are rewritten 
as 



dt 



dU_ .dU 
dx dy 



on the boundary 



(34) 



Note that we do not have enough boundary conditions to 
determine U{x, y) uniquely. In fact we can allow in Eq. 
( p^ ) arbitrary transformations of the form 



<y9 + iCz + 7 



t/i ^ -0 + 7 , tp = 7]' 



(35) 



(36) 



where C is a real constant and 7 and 7 are complex 
constants. This provides five degrees of freedom in the 
definition of the Airy potential. Two of these freedoms 
are removed by choosing the gauge in Eq.(|3^) according 
to 

dU dU 

— h i^r- = , on the boundary . (37) 

ox ay 

It is important to stress that whatever the choice of the 
five freedoms the values of the stress tensor are unaf- 
fected, and see ||^ for an exhaustive discussion of this 
point. Computing (|3^) in terms of ([l^) we arrive at the 
boundary condition 



(p{z) + zip'{z) + ■ip{z) = on the boundary (38) 

To proceed we represent ip{z) and ipi^) Laurent form: 

(p{z) = ipiz + ipo + ip^i/z + ip^2/z^ ^ , 

iP{z) = 'ijjiz + 1P0+ z + 'iJj-2 / + ■■■ ■ (39) 

This form is in agreement with the boundary conditions 
at infinity that disallow higher order terms in z. The 



The conformal map is identical in form and meaning 
to the one introduced above and successfully applied to 
mode III. On the other hand, at present we do not solve 
the Laplace equation, and our fundamental solution ( p^ ) 
is not the real part of an analytic function. We thus can- 
not simply solve in the mathematical plane and compose 
with the inverse of the conformal map. 

In terms of the conformal map we will write our un- 
known functions ip{z) and i^iz) as 

^(z)^^(<i>(«)"\z)) , ^(z)ee^($(")"\z)) . (41) 

Using the Laurent form (|l^) of the conformal map the 
linear term at cj ^ 00 is determined by Eqs. (pl|). We 
therefore can write 

i^iuj) = Vi-F^"-''^ + 1P0+ + ^^-2/^^ -I- . . . (.42) 

The boundary condition ( |38| ) is now read for the unit 
circle in the w plane. Denoting e = exp(i0) and 



t(e) = ^ (^-„/e" , v(e) 



E 

n=0 



(43) 



we write 



u(e)-H==L4^.'(6)^-^;(6) = /(6). (44) 
$'(")(e) 

The function / is a known function that contains all the 
coefficients that were determined so far: 



/(e) = -^.Ft'^ 



$'")(£) ^^Ft^ 



$'(")(e) 



(45) 



C. Solution by power series 

To solve the problem we need to compute the coeffi- 
cients (pn and '\\>n- To this aim we first represent 



$'(")(e) 

The function /(cr) has also an expansion of thc form 



(46) 



/(6)=^/.6^ 



(47) 
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In the discussion below we assume that the coefHcients bi 
and fi are known. In fact what is computed in our pro- 
cedure is the conformal map <I'^"''(a;). Thus to compute 
these coeffici ents we need to Fourier transform the func- 
tion <i>'^"^(e)/<I>'(")(e). This is the most expensive step 
in our solution, since the branch cuts that exist in Eq. 
( p^ ) rule out the use of Fast Fourier Transforms. One 
needs to carefully evaluate the Fourier integrals between 
the branch cuts. The technique how to track the posi- 
tion of the branch cuts on the unit circle was developed 
in and [? ]; after having the branch cuts the inte- 
grals are evaluated over 1000 equi-distant points between 
each pair of branch cuts. Using the last two equations 
together with (|3|) and we get 

oo 

'P-rn - ^ &-m-fc-l<^-fe = /-m , m=l,2---(48) 
fc=l 
oc 

- bm-k-i'plk = /m , m = 0, 1, 2 • • • (49) 

k=l 

These sets of linear equations are well posed. The coef- 
ficients (f-m can be calculated from equation ( ^8|) alone, 
and then they can be used to determine the coefficients 
'tjj-m- This is in fact a proof that Eq. (Q) determines the 
functions u and v together. This fact had been proven 
with some generality in 

For cracks with simple geometry this is all that we 
need. For example for a circular crack (a problem 
that was explicitly solved in ^) we simply substitute 
<I)(")((jj) — <I>(°)(cj) — tj, and proceed to solve for cp and 
■0, finding finally 



ipitu) — (fiiUJ 



(50) 



For developing cracks of arbitrary shape this is just the 
starting point. As before in the solution of mode III we 
need to compute att from which we construct the proba- 
bility measure for the first fracture event. The develop- 
ment of the then follows the same lines as before. 

To compute att at the boundary of the crack we use 
the fact that follows directly from the definitions that 



4:Re[(p'{z)] = 4Re 



$'(")(w) 



(51) 



Since this is the trace of the stress tensor, which is invari- 
ant under smooth coordinate transformation, it is also 
equal to (T„„ -|- att- Using the fact that (T„„ vanishes on 
the boundary we can write finally 



attie) = 4Re 



$'(»)(e) 



(52) 



This result is of some importance; it shows that to com- 
pute the component att of the stress tensor on the bound- 
ary we do not need to compute ipi^) at all. Of course, 
to know the stress tensor anywhere else in the body we 



need both functions. For the growth algorithm this is not 
necessary. We note that (p is computed from Eqs. (|4^), 
and this contains only &,„ with negative m. In order to 
derive a numerical scheme to compute the tangent stress 
component att on the crack we now truncate the series 
for (fi to get an approximation 



uie) 



N 
n=l 



(53) 



We see from Eq. ( [iq ) that if we wish to compute this 
series up to an order N, we need to compute the coefii- 
cients 6 , up to j < 2N + 1 and then solve the linear 
system (^8|). Note that the approximation in Eq. ( [s^ ) 
corresponds to a truncation of the series ( |4^ ) which in 
turn corresponds to a truncation of the conformal map 
Since we are interested in the macroscopic stress 
distribution along the fracture rather than in the bumpy 
microstructure, this effect is of no harm as long as we 
choose N large enough to resolve the desired patterns. 



VI. RESULTS FOR MODE I AND II 

A. Geometry without quenched disorder 

The actual fracture patterns that we find for modes I 
and II are dramatically different from those found for 
mode III for the same velocity law. 
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In Fig. 

show the fracture patterns for the linear velocity law af- 
ter about 800 fracture events. First, mode I and II are 
very similar, except for the obvious forty-five degree tilt 
in mode II due to the tilt of the symmetry axis of the 
loading. The highly ramified structure seen in mode III 
is gone, and the resulting patterns are more akin to the 
exponential velocity law in mode III, cf. Fig. ||. The 
roughening plot h{r) is also qualitatively different from 
mode III with the same velocity law. We do not observe 
a cross over to a higher exponent, indicating that there 
is no increased roughening at large scales. Indeed, for 
these modes of fracture the stress field is found to be very 
highly peaked at the tip of the fracture pattern. More- 
over, when there appear deviations towards side branch- 
ing they are quickly corrected in later growth. To make 
this point clearer we present in Fig. |l^ the stress field at 
the boundary of the crack in the vicinity of the tip. One 
can observe that the stress component is such that the 
slight tilt of the tip will be corrected at the next growth 
event. We therefore do not expect large scale roughening 
in this mode of fracture. 



B. The effect of quenched disorder 

Lastly, we present cracks with quenched disorder. First 
we followed the growth of a crack in mode I, using the 
same strategy of Subsect. 



[VC 



In Fig. show for 
example the crack obtained with ac taken from a flat 
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FIG. 11: Upper panel: fracture pattern for mode I with the 
hnear velocity law. Lower panel: Fracture pattern for mode 
II with the linear velocity law. 
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FIG. 13: The stress field at the boundary of the crack in the 
vicinity of the tip. 
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FIG. 12: The function h{r) for mode I fracture, averaged over 
11 fracture patterns. The line indicates a slope of 0.5. 



distribution with a„ 



10. Contrary to the case of 



mode III the effect of quenched disorder on the roughen- 
ing is not impressive. The roughening exponent is still 
about 0.5 for small scales, with a failure to roughen on 
the large scales. This finding remains invariant to chang- 
ing the type of quenched disorder to a power law like 
Eq. ( p2[ ) . We also do not observe roughening on the large 
scales when we put quenched disorder, and grow deter- 
ministically at the point of highest value of att — CTc. 



FIG. 14: The fracture pattern in the case of quenched dis- 
order, with (Tc taken from a fiat distribution. The pattern 
is similar to the one in Fig. |l^, with the same roughening 
behavior. 



VII. CONCLUDING REMARKS 

We have presented a solution of the problem quasi- 
static fracture using the method of iterated conformal 
map. All modes of fracture can be treated, although 
mode III is much more straightforward since the equation 
to be solved is the Laplace equation. The bi-Laplacian 
equation that is involved in mode I and II requires heav- 
ier analysis and more cumbersome numerics. Notwith- 
standing, we believe that our fracture patterns represent 
accurate solutions of the problem with the stated laws of 
evolution. 

The geometric characteristics of mode III are different 
from those of modes I and II. The fracture pattern is very 
ramified, and if we look at the backbone, (which is what 
is observed as the boundary between the two parts of the 
broken material), we find that it is rough on all scales. On 
smaller scales the roughening exponent is about 0.5, and 
on larger scales the roughening increases, having an av- 
erage roughening exponent which depends on the length 
of the fracture pattern analyzed. The exponent 0.5 is 
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intimately related to the randomness that is introduced 
by our growth rules, the higher apparent exponents are 
due to the increased ramification on the larger scales as 
is explained in Sect 



[V 



The roughening plots may ap- 
pear to be in close agreement with some experimental 
observations, which however are not conducted as mode 
III. Experimentally one expects that modes I and II are 
more relevant, but here we do not observe the cross over 
to roughness characterized by exponents of the order of 
0.75. Quite on the opposite, it appears that the rough- 
ness saturates, leading to a globally flat fracture patterns 
on the large scales. 

This leaves us with the question of how to interpret 
the observed roughness in experiments. One possibility 
is that experiments are not quasi-static. We do not have 
much to say about this possibility. Another possibility 
is that in experiments the material has remnant stresses 
and other sources of quenched disorder. This is a possi- 
bility that we can put to test. Indeed, we find that mode 
III is very sensitive to quenched disorder, cf. Subsect. 



IV C. With power law disorder we can change the geo- 
metric characteristic of the fracture patterns altogether. 
This is not the case however with modes I and II, where 



the priority of the tip in attracting the stress field is over- 
whelming. These cracks do not appear to roughen on the 
large scales even with quenched disorder. 

In summary, we believe that the experimental obser- 
vations pose an interesting riddle whose resolution will 
need a careful assessment of the experimental conditions 
and their inclusion in the theory. It is our hope that the 
solution presented above will turn out to be a useful tool 
in achieving this goal. 
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